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Abstract 

We present a novel approach to calculating strong field ionization dynamics of multielectron 
molecular targets. Adopting a multielectron wavefunction ansatz based on field-free ah initio neu- 
tral and ionic multielectron states, a set of coupled time-dependent single-particle Schrodinger 
equations describing the neutral amplitude and continuum electron are constructed. These equa- 
tions, amenable to direct numerical solution or further analytical treatment, allow one to study 
multielectron effects during strong field ionization, recollision, and high harmonic generation. We 
apply the method to strong field ionization of CO2, and suggest the importance of intermediate 
core excitation to explain previous failure of analytical models to reproduce experimental ionization 
yields for this molecule. 
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I. INTRODUCTION 



Present theoretical tools for calculating strong field ionization of atoms and molecules 



'all into two categories, 
l| and/or ADK theory 



semianalytical theories based the Strong Field Approximation 
, often with improvements over the traditional formulation to 
incorporate molecular targets 3|, |4|, and 2) direct time-dependent numerical solution of the 
Schrodinger equation. The first category suffers from approximations necessary to allow a 
semianalytical treatment, most notably the neglect of the target-specific binding potential 
of the molecular core on the ionization, continuum, and recoUision dynamics. The second 
category has the shortcoming that full numerical treatment becomes impossible as the num- 
ber of degrees of freedom increases. Time-dependent numerical solutions of the Schrodinger 
equation including a strong laser field is only feasible for one- or two-particle systems. Steps 
have been made along the numerical route to incorporate multielectron effects into strong 
field dynamics through the use of time-dependent Hartree-Fock theory {s], multiconfigura- 
tional time-dependent Hartree-Fock 6], time- dependent configuration interaction singles [tI, 
and time- dependent density-functional theory 8j. 

In this work we address both the problems of including the binding potential consistently 
throughout the strong field dynamics as well as the problem of accounting for a major frac- 
tion of multielectron effects. In particular, motivated by recent experiments demonstrating 
effects of multiple final ionic states in high harmonic generation (HHG) [9], we focus on 
a multiple ionic channel effects in strong field ionization which is the first step in HHG. 
We consider only the electronic problem, with the nuclei held fixed and work in the length 
gauge. Our approach to strong field ionization of multielectron targets combines ab initio 
quantum chemistry multielectron wavefunctions with single particle time-dependent numer- 
ical grid solutions. We use as a basis the field-free n-electron neutral and the lowest few 
(n-l)-electron singly ionized states. Any coupling to the multiply-charged ionic states is 
neglected. The wavefunction of the n*^ continuum electron associated with each ionic state 
is represented by a 3D Cartesian numerical grid. Equations of motion describing the evo- 
lution and coupling of the basis state amplitudes and the n^^ electron wavefunction are 
derived from the multielectron Schrodinger equation and contain no adjustable parameters. 
Our method is closely related to the R-matrix theory of electron-molecular scattering 10|. 



We use the identical wave function ansatz. R-matrix theory accounts for antisymmetriza- 
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tion exactly and is applicable to time-independent problems while our formalism includes 
antisymmetrization approximately but can be applied to time-dependent problems. 

As a first ^example, we apply the method to the strong field ionization of CO2. A recent 
experiment [11] found that predictions made using MO-ADK for strong field ionization of 
CO2 failed to account for the experimental angle-resolved ionization yields. Strong field ion- 
iozation of this molecule has also been theoretically analysed in recent papers using TDDFT 
in Ref. [1^ and single-channel frozen-core approach in Ref . [isl . Following our analysis pre 



sented below, we suggest that an intermediate excitation channel not considered in Ref. [11 1 
is responsible. In this channel, first an excitation of the outer-lying electron occurs concomi- 
tant with an ionic core excitation. The excited ionic core then couples back to the ground 
state of the inner core via laser coupling followed by release of the outer-lying electron. 

II. LENGTH GAUGE THEORY FOR ONE-ELECTRON CONTINUUM 
A. Hamiltonians and States 

The (non-relativistic) Hamiltonians of the laser-free ion and neutral are 



n-l 



i=l 



1 - n-l ^ 



(1) 



j=i+i 

n-l 

|r,- — r, 



where {r}„_i are the (^-1) spatial electronic coordinates of the ion, {r}n are the n spatial 
electronic coordinates of the neutral, and V^mc(^) is the electrostatic potential of the nuclei 

Vnucirl = J2 (3) 

a \r-Ra\ 

where Za and Ra are the charges and positions of the nuclei. Note that Hartree atomic units 
{h = rrie = e = 1) are used throughout. In the length gauge and dipole approximation, the 
Hamiltonian of the full n electron system interacting with the laser field is 

n 
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Let \Nj) and \Ij) be the orthogonal n-electron eigenstates of the field-free neutral and the 
(n-l)-electron eigenstates of the field-free ion respectively 



H^'lNj) = Ef\Nj) (5) 
H'\I,) = El\I,) 



Note that |A^j) and depend on both spatial as well as spin coordinates of the electrons. 
In practice, ah initio multielectron methods provide only approximate eigenstates. The 
approximate nature of \Nj) and could be taken into account by using the expectation 
value equations 

(iV,|i7^|iV,) = Ef (6) 
{I,\H'\I,) = Ej 

instead of the eigenvalue equations Eqs. ([5]). In this case, whenever a term like H^\Nj) is 
encountered in the derivation, it must be replaced by the expansion 

H'^\N,) = J2m{mH''\N,) (7) 

i 

and likewise for the terms H^\Ij). Thus additional terms coupling the basis states \Nj) and 
\Ij) will arise that are not found in the formulation when Eqs. hold. For the present work 
it is assumed that the states are the exact neutral and ionic eigenstates and Eqs. are 
used in the following derivation. In the following only the neutral ground state |A^o) = |^) 
will be used. 

B. Antisymmetrization 

We use a wavefunction ansatz that has the form (see below for specific ansatz used) 

m)) = -A\%{t)), (8) 

where |\E'p(t)) is a non-antisymmetrized 'proxy' wavefunction ansatz that treats the n*'^ 
electron differently than the remaining (n-l) core electrons. 



(9) 



is the antisymmetrization operator that antisymmetrizes the n^^ electron with the remaining 
{n-1) electrons, and Pj„ is the permutation operator that interchanges electrons j and n. 
Note that the (n-1) core electrons are already correctly antisymmetrized due to the use 
of fully antisymmetric |A^) and \I„i) states. If exact propagation of ra-electron states were 
possible and if the proxy wavefunction \^p{t)) spanned the full multi-electron space, the 
time evolution of Eq. ([8]) would be given by 

U{t,to)\^{t,)) = U{t,to)A\^p{to)) (10) 



J 

= AU{t,to)\^p{to)), 

where U(t,to) is the evolution operator defined by 

z^U{t, to) = H^ma, to), Uito, to) = T. (11) 

Equation ( |T0|) demonstrates that, at least in the case of exact propagation, one need 
not propagate a fully antisymmetrized wavefunction. Rather, it is enough to propa- 
gate a partially symmetrized initial state and apply antisymmetrization at the final time: 
AU{t,to)\^p{to)). 

With this property of time evolution in mind, we proceed to construct a propagation 
scheme for a non-antisymmetrized proxy wavefunction ansatz 

\^p{t)) = U{t,to)\^p{to)) (12) 

where the n*'' electron is treated differently than the (n-l) core electrons. The correctly 
antisymmetrized wavefunction can then be retrieved using Eq. ([8]). Since the propagator 
construct below is only approximate, due to the use of a truncated basis of ionic states, 
the reconstructed antisymmetric wavefunction will no longer be an exact representation 
of time evolution of the initial antisymmetric wavefunction. We will return to this point 
following the definition of |\l/p(t)) below to see what our propagation scheme missed using 
this procedure. 
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C. Projectors and Wavefunction Ansatz 



We wish to construct a propagation scheme based on coupled single-particle Schrodinger 
equations. With this goal in mind, we now introduce a set of single-particle orbitals that 
arise naturally for the present problem, and the multi-electron partitioning that will be used 
below. 

Given the neutral ground state |A^) and ionic states we introduce the set of (single- 
particle) orbitals, called ionization source orbitals, defined as the overlap between the neutral 
and ionic states 

10^) = iUN) (13) 

where the integration is over the (n-l) electrons of the ion. These source orbitals are related 
to the Dyson orbitals that arise in photoionization processes Q, isl by a simple 

scaling factor, = v^l0m)- addition, it will be convenient to use the normalized 

source orbitals defined as 

as well as the amplitude rj^'- 

vm = {4>i\€)- (15) 

Using 10^) and its associated ionic states \Im) we define the multi-electron source- ion states 

\Sm) = \4>i)\Im). (16) 

We now introduce the set of projectors used below to partition the multi-electron wave- 
function: 

= \S^){Sm\ (17a) 
= \N){N\ (17b) 

mAff^m (t-J2v^] 

\ k' / \ k / 



k' / \ k 



I-V'-V'^] [I-V^'-V'^) (17c) 
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where 

\N)=Arf, (/-E^-) =^N[\N)-J2'^m\Sm)] (18) 

\ m / m 

is the (normahzed) component of the neutral ground state orthogonal to the set of source-ion 
states \Sm) used, and 

is the normalization factor of the state |A^). These projectors split the multi-electron space 
into three parts with distinct physical interpretation: the project onto the overlap be- 
tween the neutral and ionic states, projects onto the component of the neutral that is 
orthogonal to all of the ionic states, and the "P^ project onto the component of the ionic 
channels that is orthogonal to the neutral. 

The projectors defined above obey the standard relations for a mutually orthogonal set 
of projectors 

pNpN ^ jyN ^20a) 
V^V^ = SmkVi (20b) 

T^LH = ^mkVL (20c) 

V^V^ = V^V^ = (20d) 

= ViVi = (20e) 
V^V^ = vlfi^ = (20f ) 

where 6mk is the Kronecker delta. Further, Using these relations it can be shown that 

{Im\V^ = ni{Iml (21) 

where TZ^ = (1 — |0m)(<^ml) projects out (removes) the source orbital from the one-particle 
space connected to the \Im) channel. Equation fl2T]) will be used below. 

The wavefunction ansatz for the proxy wavefunction constructed in the space spanned 
by these projectors is 

\^p{t)) = b{t)\N) + J2 [(^Ut)\Sm) + \XUt))] (22) 

m 

where 

\Xm{t)) = \Xm{tWm) (23) 



and \xm(t)) is the single-particle function that represents the excited n*'* electron associated 
with the ionic channel \Im), that is, \xm{t)) contains the continuum electron wavefunction 
that we wish to calculate. By imposing the condition {Sm\Xm(t)) = {(pm\Xm(t)) = 0, which 
must be enforced in the initial condition and is maintained during the propagation through 
the use of the projection operators below, the basis states in |\E'p(t)) represent orthogonal 
spaces that can be accessed by operating with the projection operators 

V^\^p{t)) = b{t)\N) (24a) 
V^\^p{t)) = am{t)\Sm) (24b) 
Kl^pit)) = IxUmim). (24c) 

Returning to the issue of antisymmetrization discussion in the previous section, we can 
now point out the dominant interactions that are neglected using the procedure 

U{t,to)A\^pito)) ^ AU{t,to)\^,{to)) (25) 

with the ansatz define in Eq. (l22l) . First we note that by using fully antisymmetric neutral 
|A^) and ionic states \Im), correct antisymmetrization is present in the (n-1) core electrons. 
Thus the procedure in Eq. (123]) only concerns the n*'^ (i.e. continuum) electron. When using 
a truncated basis of only a few low lying | /„) states, the representation given by Eq. (!22|) only 
allows for a single electron (the n^^ electron) to be in highly excited or continuum states. 
Thus, no interactions that couple a continuum (or highly excited) state of one electron 
with a continuum state of a different electron are allowed in the present formulation. Note 
that these interactions are different than interactions of two electrons simultaneously in the 
continuum, and would appear as two-particle operators that cause transitions between two- 
electron states where, for example, a continuum state of electron j and a bound state of 
electron k simultaneously couple to a continuum states of electron k and a bound state of 
electron j. 

D. Full Propagation Equations 

Consider now the Schrodinger equation for |\I/p(t)) (where dt = d/dt) 

idt\^pit)) = H''{t)\^p{t)). (26) 
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The solution of this equation is equivalent to solving ?7(t, to)l^p(^o)) discussed above. Using 
the projection operators, the Schrodinger equation becomes 

k k 
k k 
k k 

By projecting out |iV), 15*^), and l/m), and recalling Eq. fl2T|) . a coupled set of Schrodinger 
equations for am(t), h{t) and \Xmit)) is obtained 

idthit) = {N\H^{t)\N)bit) + Y,{N\H^it)\Sk)ak{t) + 5^(iV|iJ^(t)|X,(t)) (28a) 

k k 

zdtam{t) = (S„|ff^(t)|iV)6(t) + ^(^^|ff^(t)|^,)a,(t) + 5^(^„|i/^(t)|X,(t)) (28b) 

k k 

idt\xm{t)) = ni{im\H^{t)\N)b{t) + Y,ni{ijH^{t)\Sk)ak{t) + Y,ni{i^\H^ {t)\xm) 

k k 

All the required matrix elements of H^{t) are given in the Appendix. 

The set of Eqs. (128|) . together with the matrix elements appearing in the Appendix, 
is the main result of this work. In particular, they allow for the use of coupled single- 
particle propagation methods to solve for the \Xm{t)) wavefunctions rigorously coupled to 
the multielectron states |A^) and \Im)- Furthermore, numerical propagation of Eqs. fl25]) does 
not involve non-local potentials. 



III. SPECIFIC CASES AND NUMERICAL RESULTS 
A. Singlet Molecules with Uncoupled Ionic Channels 

Equations (128|) are completely general and can be applied to any target molecule regard- 
less of symmetry or charge state. In this section we chose to consider the particular case 
of ionization from singlet molecules. Further, for simplicity in the first implementation, we 
consider uncoupled ionic channels. That is, we consider ionization to multiple final ionic 
states, but calculate ionization to each individually neglecting inter-channel couplings. 

For ionization from a singlet closed-shell neutral to a particular final ion state \Im), the 
ion can be left in either spin-up or spin-down states. Thus, with spin included, every final 
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continuum-times-ion state has two spin- related channels, |/m, T) and \Lm, I), each coupled to 
a continuum electron with opposite spin, \xmit), |) and \Xmit), T) respectively. As long as 
any spin-orbit coupling is neglected, the two spin-related continuum functions are identical 
in all respects except for the differing spin label. In this case, the proxy wavefunction takes 
to form 



vl/(t)) =b{t)\N)+ al{t)\4>i,^) + \xUt),^) 



(29) 



+ 



«iWl0^,i) + IXm(t),i) |/m,T) 



and Eqs. fl28|) reduces to 



(30a) 
(30b) 
(30c) 



where am(t) = alnit) = '^Li't)j \Xm(t)) represents the (identical) spatial part of \Xmit),l) 
and \Xmit), T), and 

\<PLit)) = \Xmit)) + aUmi) (31) 

where is the (identical) spatial part of the two spin-related source orbitals t) and 
10;^,!). (In the following, we drop the explicit spin dependence of the states when the 
quantities involved to do not dependent on the spin label.) Also appearing in Eqs. (|30|) are 



\nUt)) = [Hrn - Fit) . (r- - rfl„)]|0^(t)) 



(32) 



where 



is the single-electron field-free Hamiltonian for the n*'* electron moving in the field of the 



(33) 



m*^ ionic state. 



(34) 



fc=i 



is the electronic dipole moment of the ion. 



rt-l 



(35) 
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m I 



is the electrostatic potential of the ion core electrons. The (single particle) orbital 
defined as 

= A/> [lm[E^ - m ■ dim - HmUi) - m ■ 10^)] (36) 

is the 'transfer orbital' that couples |0^)|/m) and \Xm(t))\Im) to the |A^) component of the 
neutral, where is given by 

n-l 

\^Z) = {UY.n\N)- (37) 

k=l 

This single-particle function \4>^) represents an ionization (or excitation) process where the 
laser field acts on a bound electron, but ionizes (or excites) a different electron. We refer 
to this orbital as a 'cradle orbital' in analogy with Newton's cradle, a multi-ball pendulum 
where one ball receives a force causing a different ball to swing. The remaining term in 
Eqs. ( 130|) given by 

+ F{t)- [d^ + 2\7]^\'dl^ + 2\r]^\'{4>i\fn\4>i) 

+ 2r^^{$g\4>i) + 2^(0^10^)] } (38) 

is the energy of the \N) state in the presence of the laser field, and 

n 

d^ = -{N\Y,n\N) (39) 

k=l 

is the electronic dipole moment of the neutral. The initial condition corresponding to all 
population in the neutral state are 



6(t = 0) = Vl-2|r/™|2 (40a) 

aUt = 0) = Vm (40b) 

\Xm{t = 0)) = (40c) 

The propagation equations flHU]) coupling the continuum electron \Xm{t)) to the ground 
state amplitudes am(t) and b(t) are perhaps not so transparent at first glance. They can be 
simplified in the case of negligible depletion and distortion of the ground state, 

b{t) ^ h(t = 0)e-'<W* (41a) 

am{t) ^ am{t = 0)e-^<W*, (41b) 
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where = Eq —F{t)-d^ — a\F{t)\^ takes into account a small Stark shift of the neutral. 

In this case, Eq. (l30b) simplifies to 

idt\Xn.{t)) = ^'^{[Hm - Fit) ■ (fn - dlj]\Xm{t))} (42) 

+ Fit) .[f^<Pi) + m]}e-^'^^'^'. 

This last equation is now very close to a standard laser-dressed single-particle Schrodinger 
equation for \Xm,{t))- The only difference is that orthogonality with the neutral is maintained 
through the appearance of TZ^, and the term Tl^{—F{t) ■ + l<?m)]} ^^ts as the source 

that populates \Xm(t)). For regimes where negligible depletion is expected and where Stark 
shifts and distortions of the neutral are small, Eqn. fB2l) could be used instead of Eqs. (130|) . 
In the following calculations, we use Eqs. (15UI) throughout. 



B. Ionization of CO2 



We now apply this formalism to the strong field ionization of CO2. Recently, angle- 
resolved ionization yields have been measured Q for this molecule, where the angle is 
between the molecular axis and the polarization direction of a linearly polarized laser field. 



In Ref. Ill] it was found that the experimental angular ionization pattern for CO2 differs 
strongly from the results of molecular ADK theory (MO-ADK), a single-active electron 
quasi-static tunneling theory of molecular ionization ^. The central difference is that MO- 
ADK predicts ionization peaks at an angle of ~ 30" while the measure show strong peaks 



TABLE I: Multielectron states and energies used in the CO2 ionization calculations. Zero of energy 
was set equal to the (degenerate) ionic ground state. 



State Label Energy (eV) Hole 

\N) X^Sg -13.76 

|Jl) ^^'^g,x HOMO 

II2) X^Hg^j^ HOMO 

\h) A'^Ilu,x 3.53 HOMO-1 

IJ4) A^I^u,y 3.53 HOMO-1 

1/5) 4.28 HOMO-2 
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FIG. 1: (Color online) Angle-resolved ionization yields for intensity 1.5 x 10^"^ W/cm^. 9 is angle 
between the molecular axis and the polarization axis of the laser field. The solid curves are results 
using a grid spacing A = 0.1 a.u. while the dashed curves used A = 0.2 a.u. 



at ~ 45°. Here we consider angle-resolved ionization yields of CO2 exposed to a single cycle 
of an 800 nm laser {oj = 0.057 a.u.). 

The neutral |A^) and lowest five ionic \Im) multielectron orbitals are calculated using the 



GAMESS electronic structure code 



16[ . All calculations use the cc-pVTZ basis set 17| and 



were done at a CAS level using 16 (neutral) or 15 (cation) active electrons in 10 orbitals. 
Experimental geometry of the CO2 ground state is used (linear, i?c-o=l-1621 A). The 
states and energies used are shown in Table [T] along with the approximate location of the 
hole (relative to the neutral) left by the removed electron for each ionic state. Equations 
( |29|) are solved using a leapfrog algorithm. The wavefunction \Xm(t)) is represented on a 
3-dimensional Cartesian grid. The grid extends to ± 13 a.u. in the x and z directions, and 
to ± 8 a.u. in the y direction. All calculations are done with a grid spacing of A = 0.1 
a.u. in all directions unless otherwise specified. Absorbing boundary conditions are used 
in the xz plane with a width of 5 a.u. from the boundary edges 18|]. The ionization yield 
was calculated by monitoring the density absorbed at the boundaries. The CO2 molecule 
has the C atom at the origin and has the bond axis aligned along the z-axis. The laser 
field F{t) = Sq sin{ujt) is rotated in the xz plane. The angle 6 is the angle between the 
laser polarization and the molecular axis. The laser field is turned off after a single cycle, 
F{t > 2tt/uj) = 0, and the simulations are run until t = 150 fs (an additional 40 fs after the 
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FIG. 2: (Color online) Left column: Dyson orbitals of the ionic channels considered. Angle between 
the laser field and the molecular axis is depicted in the top left panel. Only the Dyson orbitals 
with lobes in the plane of the laser field are shown. Center and right columns: Angular ionization 
yields for intensity 1.5 x 10^^ W/cm^ for each ionic channel considered. Bottom-right panel shows 
the total ionization yield which is the sum of all channels. 

single cycle is over) to allow the liberated electron density to be absorbed at the boundary. 
A time step of At = 0.00133 a.u. is used for the time propagation. 

Figured] plots angle-resolved ionization yields for the five final ion states considered for a 
intensity of 1.5 x 10^^ W/cm^. The solid lines correspond to calculations with the step sizes 
specified above, while the dashed lines show results using A = 0.2 and At = 0.00266 a.u. 
While the total yields continue to decrease a bit as the grid size becomes finer, the general 
character and relative behavior of the ionization channels is preserved. For all angles, the 
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ionization yield is dominated by the X^IIg channels. Polar plots showing the angular shape of 
each ionization channel are presented in Fig. [2l Also shown in this plot is the total ionization 
yield that included the yield from all channels (bottom-right panel), which is effectively the 
same as the yield including only the two X^IIg channels (not shown). The total ionization 
yield has a 'bow tie'-like pattern, with peak values appearing near 30°. This is in closer 
agreement with the MO-ADK results than the experimental distributions, both presented 
in Ref. (llj. Note that the MO-ADK results of Ref. (llj include only the 'in-plane' HOMO 



channel which would correspond to the X^Ug^x channel alone. Thus, our uncoupled channel 
calculations still fail to reproduce the experimental peak positions seen in Ref. 



C. Role of Nodal Planes and the Binding Potential 

It has been shown that the presence of nodal planes in the ionizing orbitals leads to 
suppression of the ionization rate {sl. Most prominently, large suppression is expected to 
occur when the laser field is aligned along a nodal plane. This expected trend can be seen 
in our results (Fig. [2]) by comparing the angular ionization yields with the corresponding 
Dyson orbitals. However, two features stand out that deserve attention. First, although 
suppression is seen along both nodal planes in the X^H^ ^.^ distribution, there is much more 
suppression along the 90° node than along the 0° node. Second, there is a dip in the A^H^^^; 
ionization yield at 90° that corresponds to no obvious feature in the A'^Uu^x Dyson orbital. 

Consider the X^H^, 2, distribution. We first consider the case when the peak laser field is 
1 X 10^^ W/cm^ and return to the case of 1.5 x 10^^ W/cm^ below. (Note that although 
the angular ionization yields shown in Fig. [2] where calculated for 1.5 x 10^^ W/cm^, the 
angular shapes for each channel are very similar when using an intensity of 1 x 10^^ W/cm^.) 
Panels (a) and (b) in Fig. [3] plot the X^Hg^^,' Dyson orbital along with select contours of the 
instantaneous potential at the peak of the laser pulse for an intensity of 1 x 10^^ W/cm^. 
The contours are taken at the ground state energy of the neutral and show the entrance and 
exit of the tunneling barrier through which the X^Hg .j. Dyson orbital must escape. Panel 
(a) shows the contours when the laser is aligned along 0°, while panel (b) is for the 90° case. 
The short solid lines connecting the entrance and exit depict the tunneling path positioned 
along the peak of the orbitals lobes. The tunneling barrier along these paths are shown in 
panel (c). Already one can see that the tunneling path through which the orbital lobes must 
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FIG. 3: (Color online) Top row: Dyson orbital of the ^Tig^x channel and ground state energy 
contours of the potential energy landscape at the peak of the laser field for 1 x 10^^ W/cm^. In (a) 
the laser points along 0° while in (b) the laser field points along 90°. Also shown are the tunneling 
paths through which the orbitals lobes must tunnel. Panel (c) plots the potential energy along 
the tunneling paths. The solid line denotes the ground state energy of -13.76 eV = -0.5058 a.u. 
Panels (d) and (e) plot the X^IIg^^,' Dyson orbital along with the ground state energy contours of 
the potential energy for 1.5 x 10^^ W/cm^. In panel (d) the laser field points along 0° while in 
panel (e) the laser field points along 90". 



pass is much shorter in the 0° configuration than in the 90° configuration suggesting the 



origin of the difference in suppression in the 0° and 90° degrees directions seen in the X^II 



distribution. In order to get a quantitative semiclassical estimate of the ratio of ionization 
at 0° and 90°, we use the WKB tunneling formula 



Rate ~ exp 



' v/2(V(x') - E)dx' 



XO 



(43) 



16 



10 
5 



m 



-5 
-10 





(a) 


*' ^ ■ ■ ■ 



-10 



-5 5 
z (au) 




10 

10 -10 



FIG. 4: (Color online) Dyson orbital of the A'^Hu^x channel and ground state energy contours of 
the corresponding potential energy landscape at the peak of the laser field for 1.5 x 10^^ W/cm^. In 
(a) the laser points along 90° while in (b) the laser field points along 45°. The thick line segments 
lines show the shortest the tunneling paths. 

where the integral is taken across the tunneling barrier. Using this measure, we find that 
the rate of tunneling along 0° should be larger than the rate along 90° by a factor of 7.3, 
which in good agreement with the actual ratio of 8.8 extracted from the simulations. Panels 
(d) and (e) plot the same contours as in panels (a) and (b), but now for the intensity 
of 1.5 X 10^^ W/cm^. In this case, the ionization is above barrier, and the ground state 
energy contours show the 'doorway' opened by the presence of the laser field. Although 
a quantitative estimate is difficult in the above-barrier regime, one can see that for 0° the 
doorway encompasses almost the whole width of the Dyson orbital along this direction, while 
for 90° the doorway is allowing only a small portion of the orbital localized around the nodal 
plane to pass. Thus, the analysis of the potential landscape in the 1.5 x 10^^ W/cm^ case 
allows for a qualitative understanding of the large difference in suppression along the nodal 
planes seen in the X.'^Hg^x angular ionization yields. 

We turn now to the A^n^ .j, channel, where a similar analysis accounts for the dip at 90°. 
Figure H] shows the A'^Uu^x Dyson orbital along with the ground state energy contours. In 
panel (a) the laser field points along 90° while panel (b) corresponds to 45°. Both panels 
correspond to a peak intensity of 1.5 x 10^^ W/cm^, the case shown in Fig. [21 Integrating the 
tunneling rate along the paths shown in the plots, which are the shortest paths connecting 
the inner and outer regions in both cases, we calculated that the tunneling rate for the 90° 
case should be suppressed by a factor of 0.6 as compared to the 45° case. This is again in 
good agreement with the actual suppression of 0.5 extracted from the results in Fig. [2] for 
this channel. 
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D. Toward Coupled-Channel Ionization of CO2 



We can use the results of the present uncoupled channel calculations to infer potentially 
important ionization mechanisms that will appear in a coupled channel treatment. 

In our formulation, the wavefunction \Xm{t)) carries not only the continuum states, but 
also a complete set of bound states bound to the \Im) ionic core. Thus, using the same 
simulations discussed above, we can calculate excited, but un-ionized, population of n*'^ 
electron surrounding each ionic core. In particular, the top two panels of Fig. [5] show the 
angular excitation yields surrounding the A^n„ ,j, ionic core for two intensities of IxlO^"^ 
and 2.5x10^^ W/cm^. These yields show strong peaks near (or beyond) 45°. In addition, 
as shown in Fig. [H^a), the peak excitation yield surrounding the A'^U^ .j. ionic core is much 
larger than the peak ionization yield coming from the X^II^ channels. 

In an uncoupled channel formulation, as is the case with the present calculations, this 
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FIG. 5: Top panels: Excitation yield surrounding the A^H^j^^^ ionic core for the two intensities of 
1 X 10^^ W/cm^ and 2.5 x 10^^ W/cm^. Bottom panels: Ionization estimates for the two intensities 
showing the direct ionization channel (thick-dashed) analogous to the " sum over all channels" panel 
shown in Fig. 2, the estimated intermediate excitation channel (thin), and the sum of these two 
channels (thick). 
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FIG. 6: Panel (a): Maximum of A'^U^^x excitation and maximum of X^IIg ionization as intensity 
is varied. Panels (b) and (c) plot the population of the X^IIg ionic state using the 2-level laser 
coupled model of Eq. (j44p . with all population initially in the A'^ILu^x state. Panel (b) is for a single 
cycle pulse, and panel (c) is for a longer pulse with smooth Gaussian envelope. 

excited population surrounding the A'^Uu^x ionic core is trapped. (We have checked that 
similar excited population exists at the end of a 5 fs Gaussian laser pulse in addition to 
the single cycle pulses used herein). However, in a coupled channel formulation, some of 
this excited population surrounding the A'^Uu^x core will be moved to the X^n^^^. ionic 
core through laser-induced dipole coupling of the A'^U^^x and X^IIg .j. core, i.e. through the 
polarization of the ionic cores. The amount of ionic core coupling can be estimated by 
solving a 2-state problem for the laser coupling of the A^II^ .j. and X^IIg .j. cores 
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Cxit) 
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Cxit) 

CA{t) 



(44) 



Ea -F{t)fXAB 
-F{t)fiAB ^ ^ 

where [iab = -0.46722 a.u. is the transition dipole between the ionic states y^Iig^x and 
A^n^a;, calculated using GAMESS as outlined above, and Cxit) and C^(t) are amplitudes 
of the X^IIg .J. and A^n„ ,j, states. Figures [6] (b) and (c) plot |Cx(t)| as a functions of time for 
two different pulse, panel (b) uses a single cycle and panel (c) uses a multi-cycle pulse with 
Gaussian envelope, with the initial condition Cx{t) = and CA{t) = 1. The calculations 
were done for two different intensities, 1 x 10^^ W/cm^ (thick lines) and 2.5 x 10^^ W/cm^ 
(thin lines). These calculations allows us to estimate that about 5 to 10% of the excited 
population surrounding the A^H^ core will couple back to the X^Ug^x state on subsequent 



19 



cycles. Some (and perhaps all) of this excited population will escape the core region once 
coupled back to the X'^Ug^^ ionic core. We thus anticipate two important ionization channels 
in a coupled-channel formulation of CO2, the direct channel and an intermediate excitation 
channel 

Direct : C02{X^^g) ^ CO+iX'^Eg) + e" (45) 

Inter. Ex. : C02{X^^g) ^ C0^(A2S„,^)(e-)* COUX'^Jlg) + e~ (46) 

where (e~)* denotes an excited electron. Assuming that all of the excited population will 
escape the core upon coupling from the A^fl^ .j, back to the X^Ug .^ state, the intermediate 
excitation channel will carry predominantly the angular imprint of the C02(X^Sg) — > 
C0^(A^Su^2:)(e")* excitation step. 

The direct channel yield and (estimated) intermediate channel yield, as well as their 
sum, is plotted in the bottom two panels of Fig. [5] for the two intensities shown. Here the 
intermediate excitation channel yield was estimated by multiplying the yields for excitation 
on the A^n^,^. ionic core by the amount of coupling seen in Figs. [6](b) and (c), 0.05 in 
the case of 1 x 10^^ W/cm^ and 0.1 for 2.5 x 10^^ W/cm^. At the higher intensity, the 
direct channel dominates, while for the lower intensity the intermediate excitation channel 
is becoming important. Further, the peak of the total ionization estimate for 1 x 10^^ 



W/cm^ is now approaching 45°, as seen in the experiment llj. Our treatment of the 



proposed intermediate excitation channel is admittedly crude, and fails to reproduce the 



sharpness of the experimental peaks seen in Ref.[ll|. An accurate description requires a full 
coupled-channel treatment of Eqs. fl25]) that includes at least the X'^Ug^ and A'^U^^x states. 
However, from the scaling of the excitation and ionization yields seen in Fig. El^a) it is clear 
that the intermediate excitation channel will become important for a correct description of 
strong field ionization of CO2 at intensities up to (and perhaps beyond) 10^^ W/cm^. 



IV. SUMMARY 



In this work we developed a method for strong field one-electron ionization of multielec- 
tron targets. Our method uses field-free multielectron orbitals to describe the neutral and 
lowest few ionic states. These multielectron basis states are supplemented with a one-particle 
numerical grid used to represent the continuum electron. Equations of motion coupling the 
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basis states to the continuum grid are derived from the multielectron Schrodinger equation. 
The result is a coupled set of single-particle Schrodinger equations describing ionization into 
each final ion state included in the ionic basis. Our equations are general and applicable to 
strong field ionization of any small molecule. 

As an example, we studied ionization of CO2 in the uncoupled channel approximation 
including the lowest five ionic states of CO^. Strong field ionization of this molecule has 
been experimentally shown to deviate from the predictions of MO-ADK, a single-active- 
electron quasi-static model of molecular ionization. Our method allows the inclusion of 
two dominant effects not present in MO-ADK: 1) influence of the specific shape of the 
tunneling barrier discussed in Sec.III-C and 2) the possibility to rigorously couple multiple 
ionic channels as dissused (but presently not implementd) in Sec.III-D. In our analysis, 
the deviations from MO-ADK seen experimentally likely arise from intermediate ionic core 
excitations followed by interchannel coupling. 



APPENDIX A: MATRIX ELEMENTS OF THE HAMILTONIAN 



In order to evaluate the matrix elements appearing in Eq. fl28p . we need to know how 
H^{{fn},t) acts on the basis states. The Hamiltonian acting on the neutral state gives 



H^\N) = {E^-J2F{t)-rA\N). 



(Al) 



fc=i 



The Hamiltonian acting on a state \(f)m)\Im), where \Im) is an ionic state and is here 
an arbitrary single particle function, gives 

n-l 

I • — , ^ 

rk \(j)m)\Im) 



{\(j)m)\Im)) 



^ n— 1 ^ n 



k=l ' k=l 

n-l 



:fk — f„ 
k = l ' " k = l 



rk 



(A2) 



\Im) 



k=l 



where 



{H^fMm)) \Im) + I E VfjfnMm) ] ' ^(t) ■ ( J^i^J.m ' 4j\<Pm) ] |/,) 



(A3) 
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is the single-electron field-free Hamiltonian for the n^^ electron coupled to the ionic state 

\Im) t 

n-1 
k=l 

are the electronic dipole moments and transition dipoles of the ionic states, and 



n— 1 ^ 
k=l 

are the electrostatic potentials and inter-ionic couplings. Eq. (IA2I) is only exact if a complete 
basis of \Ij) is used. If this basis is truncated, Eq. (1A2I1 gives {\4>m)\Im)) projected into 
the space of the truncated basis. Below we will also need the electronic dipole of the neutral 
defined as 

n 

= -{N\Y,n\N) (A6) 

k=l 

We now calculate the required matrix elements of the Hamiltonian. First First consider 
the matrix elements of the 'primitive' basis functions |A^), \Sm), and \Xm{t)). In the following 
matrix elements the convention m 7^ is used in order to avoid excessive use of Kronecker's 
delta. 

{Srn\H^{t)\Sr^) = {{4>i\{Irn\)H^m\4>i)\U) 

= {4>i\Hm\4>i) - m ■ [{&n\~^i) - 4J (A7) 

{S^\H^{t)\S,) = {{^i\{Irn\)H^m~^l)\h)) 

= {^XSD + Fit) . d^^S'M) (A8) 
{SjH^{t)\N) = {{4>i\{Im\)H^{t)\N) 

= VmE^ - Fit) . [Vmi&nWl) + i^HfZ)] (A9) 
{Sm\H^it)\XUt)) = i{4>i\{Im\)H^it)i\xUtWm)) 

- '"^^ \Hm\Xn.it)) - Fit) ■ (0^|r-;|x™(t)) (AlO) 



m I 



{S^\H^it)\X,it)) = i{4>i\{Im\)H^{t)i\Xkit))\h)) 

= {4>X,\Xk{t)) + F{t) ■ dl,{4>i\xk{t)) (All) 
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{N\H^{t)\N) = E^ + F{t)-d^ 
{N\H^{t)\Xr,{t)) = {NlH^mXmitWm)) 

= -F{t) ■ [v*m{4>i\rn\Xm{t)) + {^IxUt))] 



(A12) 



(A13) 



{UH^{t)\s^) = {UH^m<pi)\ij) 

= Hm\4>i)-F{t)-{fn-dlMi) 



(AM) 



(A15) 



{I^\H'{t)\N) = VmE^\<P'j - F{t) ■ [fn\<P'jVm + 

{i^\H^{t)\x^{t)) = {i^iH^mxmmim)) 

= Hm\Xn.{t)) - F{t) ■ {rn - d'^JWmit)) 



(A16) 



(A17) 



(A18) 



= V^,\xk{t)) + F{t)-dl,\xk{t)) 

Now we use these matrix elements to evaluate the remaining terms in Eqs. (!28ll that 
involve |A^) 



{S^^\H^{t)\N) = AfJ{S:^\H^{t)\N) -J2vk{Sm\H^{t)\Sk) 



(A19) 



N 



{S^\H^{t)\N) - r]m{Sm\H''{t)\Sm) - J2 Vk{Sm\H^{t)\Sk) 



N 



VmE^ - Fit) ■ [Vmi&ml^l^'J + 



S \^ |I5 ' 



IS \TC\ 



s iTC\ 



-? 1 

"mmj 



ky^m 
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{N\H^(t)\Sm) = {Sm\H^(t)\Nr 

= [v*jE^ - i^ilHmm] - F{t) ■ [&D + v*md: 



(A20) 



^ 1 

mm J 



\ m / \ k 

mk 

+ \^fN\'J2 [l^mr(-5„|i^^WI-S„) -r7m(A^|i^^(t)|-Sm) -r7;(^m|i^''W[i^^jl) 



{N\H^{t)\X^{t)) = AfJ{N\H^{t)\Xrr.{t))-J2vl{Sk\H^{t)\Xr.{t)) 



N 



- Fit) ■ KA&n\Xmit)) + i&mm 



- v*J{<p'jHMUt)) - F{t) ■ (0^|r-;|xmW)] 

- E vi[{4>x\xkit)) + m ■ 4nk{4>i\xkm 

k 

- E vl[{4>i\V^k\Xk{t)) + Fit) . diSilXkim] (A22) 



{Im\H^it)\N) = A/-^[(7^|//^(t)|Ar)-^r7fe(/„|//^(i)|5fe) 



^o''|0^>^m - i^(^) • rn\4>i)Vm - Fit) ■ |0 

~^)+r]^Fit)-ifn-dlj\^i) 



kj^m 



AA^ [rim[E^ - Fit) ■ diJl^i) - VmHml^i) - Fit) 

kj^m 



(A23) 



24 



[1] L.V. Keldysh. Zh. Eksp. Teor. Fiz. 47, 1945 (1964) [Sov. Phys. JETP 20, 1307 (1965)]; F.H.M. 

Faisal. J. Phys. B 6, L89, (1973); H.R. Reiss. Phys. Rev. A 22, 1786 (1980). 
[2] N.B. Delone and V.P. Krainov, Multiphoton Processes In Atoms, 2'*'^ ed., Springer-verlag 

(1994). 

[3] J. Muth-Bohm, A. Becker, and F. H. M. Faisal, Phys. Rev. Lett. 85, 2280 (2000). 
[4] X.M. Tong, Z.X. Zhao, and CD. Lin, Phys. Rev. A 66 033402, (2002). 

[5] K.C. Kulandcr, Phys. Rev. A 36, 2726 (1987); M. Lein, N. Hay, R. Velotta, J.P. Marangos, 

and P.L. Knight, Phys. Rev. Lett. 88, 183903 (2002). 
[6] M. Kitzler, J. Zanghellini, Ch. Jungreuthmayer, M. Smits, A. Scrinzi, and T. Brabec, Phys. 

Rev. A 70, 041401(R) (2004); T. Kato and H. Kono, J. Chem. Phys. 128, 184102 (2008). 
[7] N. Rohringer, A. Gordon, and R. Santra, Phys. Rev. A 74, 043420 (2006). 
[8] D.A. Telnov and Shih-I Chu, Phys. Rev. A 79, 041401(R) (2009). 

[9] W. Li, X. Zhou, R. Lock, S. Patchkovskii, A. Stolow, H.C. Kapteyn, M.M. Murnane, Science 
322, 1207 (2008); B.K. McFarland, J.P. Farrell, P.H. Bucksbaum, and M. Guhr, Science 
322, 1232 (2008); O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. 
Corkum, M.Yu. Ivanov, Nature 460, 972 (2009). 

[10] See, e.g., P.O. Burke and J. Tennyson, Molecular Physics 103, 2537 (2005). 

[11] D. Pavicic, K.F. Lee, D.M. Rayner, P.B. Corkum, D.M. Villeneuve, Phys. Rev. Lett. 98, 
243001 (2007). 

[12] Sang-Kil Son and Shih-I Chu, Phys. Rev. A 80, 011403(R) (2009). 
[13] M. Abu-samha and L. B. Madsen, Phys. Rev. A 80, 023401 (2009). 
[14] B.T. Pickup, Chem. Phys. 19, 192 (1977). 
[15] Y. Ohrn and G. Born, Adv. Quantum Chem. 133, 1 (1981). 

[16] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, 
N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, and J.A. Montgomery, J. 
Comput. Chem. 14, 1347 (1993). 

[17] T.H. Dunning, Jr., Chem. Phys. 90, 1007 (1989). 

[18] D.E. Manolopoulos, J. Chem. Phys. 117, 9552 (2002). 



25 



